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This  final  report  summarizes  our  major  accomplishments  on  the  research  under  AFOSR  Grant  FA9550- 
13-1-0092  during  the  project  period  between  March,  2013  -  March,  2016. 

In  this  project,  we  continue  our  development  of  our  Riemann-solver-free  spacetime  discontinuous  Galerkin 
method  for  general  conservation  laws  to  solve  compressible  magnetohydrodynamics  (MHD)  equations.  The 
method  is  first  applied  to  solve  the  3x3  MHD  model  system  in  the  phase  space  which  exactly  preserves  the 
MHD  hyperbolic  singularities.  Numerical  results  show  that  the  method  is  able  to  solve  the  model  system 
correctly,  which  makes  the  method  very  promising  in  solving  the  complete  ideal  MHD  equations.  The  method 
is  then  extended  to  solve  the  1-D  7x7  and  2-D  8x8  MHD  equations.  The  Powell’s  approach  by  adding 
appropriate  source  terms  is  adopted  to  handle  the  V  •  B  =  0  condition.  Again,  the  numerical  results  show 
that  the  present  method  is  able  to  resolve  the  complex  MHD  waves  without  the  need  of  any  type  of  Riemann 
solvers  or  other  flux  functions.  The  success  of  solving  MHD  equations  further  strengthens  our  belief  that 
the  DG-CVS  is  an  effective  approach  in  solving  systems  where  accurate  and  reliable  Riemann  solvers  are 
difficult  to  design. 

1  Background 

This  project  continues  our  previous  AFOSR  projects  (Grant  No.  FA9550-08-1-0122  and  Grant  No.  FA9550- 
10-1-0045)  to  extend  and  verify  our  high  order  spacetime  discontinuous  Galerkin  cell  vertex  scheme  (DG- 
CVS)  toward  solving  the  magnetohydrodynamics  (MHD)  equations. 

The  magnetohydrodynamics  (MHD)  equations  are  often  used  to  model  electrically  conducting  fluid  flows 
where  the  electromagnetic  forces  can  be  of  the  same  order  as,  or  even  greater  than,  the  hydrodynamic  forces. 
Plasmas  in  hypersonic  and  astrophysical  flows  are  one  of  the  most  typical  examples  of  such  conductive  fluids. 
Though  MHD  models  are  a  low  order  approximation  of  actual  plasma  processes,  they  have  been  successfully 
applied  to  simulate  many  important  plasma  processes.  Today,  the  MHD  models  remain  powerful  tools  in 
helping  researchers  to  understand  the  complex  physical  processes  in  the  geospace  environment. 

For  example,  the  ideal  MHD  equations  combines  the  Euler  equations  of  gas  dynamics  with  Maxwell  equa¬ 
tions  of  electromagnetics  for  problems  in  which  viscous,  resistive  and  relativistic  effects  can  be  neglected.  The 
ideal  MHD  equations  contain  a  set  of  nonlinear  hyperbolic  equations  as  shown  in  the  following  conservative 
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form: 
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Here,  p,  v,  E,  P,  po  and  B  represent  the  mass  density,  the  velocity  vector,  the  total  energy,  the  hydrodynamic 
pressure,  permeability  of  vacuum  and  the  magnetic  field.  In  addition  to  satisfying  the  induction  equation 
(Eq.  (Id)  above),  the  magnetic  field  is  also  divergence  free,  i.e.  V  •  B  =  0  which  can  be  considered  as  an 
extra  constraint. 

As  can  be  seen  in  Eqs.  (la-ld),  the  MHD  equations  are  complicated  nonlinear  partial  differential  equa¬ 
tions  (PDEs)  that  require  some  numerical  method  to  solve  for  general  problems.  Since  MHD  flow  fields  may 
develop  shock  waves,  contact  discontinuities  and  shear  layers,  the  numerical  schemes  are  required  to  be  of 
high  resolution  which  are  able  to  resolve  all  kinds  of  discontinuities  accurately  and  robustly. 


1.1  Traditional  Numerical  Methods  for  MHD  Equations 

Due  to  its  close  mathematical  similarity  to  the  Euler  equations  of  gas  dynamics,  the  MHD  equation  system  is 
often  solved  by  various  methods  extended  from  those  originally  designed  for  Euler  equations.  The  Godunov- 
typed  methods  use  exact  or  approximate  Riemann  solvers  to  provide  numerical  inter-cell  spatial  fluxes.  Since 
the  MHD  equations  contain  much  more  complicated  wave  structures  than  the  compressible  Euler  equations, 
the  (approximate)  Riemann  solvers  for  Euler  equations  must  be  sophisticatedly  modified  or  even  redesigned 
to  handle  the  degeneracies  and  instabilities  of  the  MHD  equations. 

Two  of  the  most  widely  used  approximate  Riemann  solver  families,  the  Roe  scheme  [  ]  and  the  HLL 
(Harten,  Lax  and  van  Leer)  scheme  [2]  and  its  two-state  variant,  HLLC  scheme  [3,  4],  have  gained  increasing 
popularity  in  solving  the  MHD  equations.  For  example,  see  [5,  6,  7]  for  MHD  solvers  based  on  the  Roe’s 
scheme  and  [8,  9,  10,  1  ]  for  MHD  solvers  based  on  the  HLL  or  HLLC  scheme.  Roe’s  scheme  requires  eigen- 
decomposition  and  becomes  very  complicated  in  MHD  equations.  For  example,  the  Roe  averages  of  quantities 
are  not  clearly  defined  in  MHD  system.  Moreover,  due  to  the  complexity  and  non-strict  hyperbolicity  of  the 
MHD  system,  the  validity  of  the  eigensystems  of  the  MHD  system  is  not  unanimously  agreed  upon  among 
researchers.  The  HLL  or  HLLC  based  Riemann  solvers  exhibits  several  advantages  over  Roe’s  scheme.  For 
example,  the  HLLC  scheme  preserves  positivity,  satisfies  the  entropy  condition  and  does  not  need  eigen- 
decomposition  of  the  system.  The  accuracy  of  MHD  solvers  is  dependent  on  the  underlying  Riemann  solvers 
[  .2].  Therefore,  the  use  of  Riemann  solvers  adds  uncertainties  to  the  accuracy  of  MHD  solvers.  The  kinetic 
MHD  scheme  [13,  14,  15]  is  a  flux  vector  splitting  method  with  the  consideration  of  particle  transport.  It 
is  claimed  to  be  more  efficient  than  other  Riemann- solver  based  MHD  schemes  due  to  the  simplicity  of  the 
kinetic  flux  functions  [14]. 

In  addition  to  upwind  schemes,  the  fully  discrete  [16]  or  semi-discrete  [V  ]  central  schemes  can  also  be  used 
for  MHD  equations.  Fully  discrete  central  schemes  use  staggered  grids  to  evolve  cell  averages  and  eliminate 
the  need  for  the  information  of  the  eigenstructure  of  the  system.  Moreover,  central  schemes  eliminate  the 
need  for  dimensional  splitting  [16].  Central  schemes  are  considered  as  Riemann- solver- free  methods.  Central 
schemes  are  very  attractive  in  that  the  notorious  pathological  phenomena  associated  with  some  Riemann 
solvers  are  absent  in  Riemann-solver-free  methods. 

Another  issue  when  solving  the  MHD  equations  using  numerical  methods  is  how  to  preserve  V  •  B  =  0 
for  the  magnetic  field.  Failure  to  enforcing  this  divergence- free  constraint  may  lead  to  a  nonlinear  numerical 
instability.  Several  approaches  are  available  to  preserve  this  constraint.  One  approach  is  to  use  the  projection 
method  [18]  by  solving  a  Poisson  equation  for  the  magnetic  potential  which  can  be  used  to  correct  the  non¬ 
divergence-free  magnetic  field.  In  Ref.  [19],  a  DG  method  using  locally  divergence-free  piecewise  polynomials 
as  the  solution  space  is  presented  to  solve  the  MHD  equations.  Another  so-called  Powell  approach  [7]  is  to 
add  source  terms  proportional  to  V  •  B  to  the  ideal  MHD  equations.  Despite  some  success,  this  approach 
leads  to  a  non-conservative  formulation  which  may  result  in  incorrect  weak  solution.  It  is  also  possible  to 
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use  a  generalized  Lagrange  multiplier  method  [20]  to  “clean”  the  divergence  error.  Instead  of  solving  an 
elliptic  Poisson  equation,  this  approach  solves  a  damped  hyperbolic  equation  for  the  divergence  error  with 
appropriately  chosen  parameters.  In  Ref.  [21],  the  magnetic  vector  potential  A  is  solved  in  place  of  the 
magnetic  field  B  itself.  Since  V  •  (V  x  A)  =  0  is  guaranteed  mathematically,  solving  A  automatically  leads 
to  V  •  B  =  0.  However,  setting  the  boundary  conditions  for  the  vector  potential  is  not  straightforward. 
Another  very  popular  method  is  the  constrained  transport  (CT)  method  of  Evans  and  Hawley  [22].  The  CT 
method  can  be  considered  as  a  predictor-corrector  approach  for  the  magnetic  field  [23] .  It  uses  the  computed 
conservative  variables  to  approximate  the  electric  field  which  is  used  to  update  the  magnetic  vector  potential. 
The  magnetic  vector  potential  is  then  used  to  compute  the  divergence-free  magnetic  field.  The  CT  method 
does  not  need  to  solve  a  global  elliptic  equation  and  has  no  free  parameters.  The  CT  method  and  its  variants 
have  been  investigated  by  many  researchers [24,  25,  23]. 

In  this  project,  we  solve  the  following  system  [26,  7]. 
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1.2  Riemann-Solver-Free  DG-CVS  for  General  Conservation  Laws 

A  novel  high  order  discontinuous  Galerkin  cell-vertex  scheme  (DG-CVS)  [27,  28,  29,  30,  31,  32,  33]  has  been 
developed  for  hyperbolic  conservation  laws  under  the  support  from  AFOSR.  DG-CVS  was  inspired  by  the 
spacetime  Conservation  Element /Solution  Element  (CE/SE)  [34]  method  and  the  discontinuous  Galerkin 
(DG)  [31]  method.  The  main  features  of  DG-CVS  are  summarized  as  follows: 

•  based  on  space-time  formulation.  Therefore,  DG-CVS  automatically  satisfies  the  Geometric  Conserva¬ 
tion  Law  (GCL)  and  solves  moving  mesh  problems  [36]  without  special  care  as  needed  by  ALE-based 
methods. 

•  arbitrary  high- order  accuracy  in  both  space  and  time.  Space  and  time  are  handled  in  a  unified  way 
based  on  space-time  flux  conservation  and  high-order  space-time  discontinuous  basis  functions.  This  is 
in  contrast  to  semi-discrete  methods  where  the  temporal  order  of  accuracy  is  limited  by  the  Backward 
Difference  Formula  (BDF)  or  the  multi-step  Runge-Kutta  method. 

•  Riemann- solver  free.  DG-CVS  does  not  need  a  (approximate)  Riemann  solver  to  provide  numerical 
fluxes  as  needed  in  finite  volume  or  traditional  DG  methods.  The  Riemann-solver-free  feature  offers 
two- fold  advantages.  First,  this  Riemann-solver-free  approach  eliminates  some  pathological  behaviors 
(e.g.,  carbuncle  phenomenon,  expansion  shocks,  etc.)  associated  with  some  Riemann  solvers.  Second, 
it  is  suitable  for  any  hyperbolic  PDE  systems  whose  eigenstructures  are  not  explicitly  known.  The 
DG-CVS  based  solvers  have  been  successfully  applied  to  solve  scalar  advection-diffusion  equations  [31], 
compressible  Euler  equations  [27,  33],  shallow  water  equations  [37]  and  the  level  set  equation  [38]  (see 
Fig.  ??  for  simulation  examples  using  the  DG-CVS  method). 

•  reconstruction  free.  DG-CVS  solves  for  the  solution  and  its  all  spatial  and  temporal  derivatives  simul¬ 
taneously  at  each  space-time  node,  thus  eliminating  the  need  of  reconstruction. 

•  suitable  for  arbitrary  spatial  meshes.  The  CE  and  SE  definitions  in  DG-CVS  are  independent  of  the 
underlying  spatial  mesh  and  the  same  definitions  can  be  easily  extended  from  1-D  (cf.  Fig.  1)  to 
higher-dimensions  (cf.  Fig.  2)  without  any  ambiguity.  Though  not  shown,  the  same  CE  and  SE 
definitions  also  apply  to  meshes  with  hanging  nodes  [30]. 

•  highly  compact  regardless  of  order  of  accuracy.  Only  information  at  the  immediate  neighboring  nodes 
will  be  needed  to  update  the  solution  at  the  new  time  level.  Compactness  eases  the  parallelization  of 
the  flow  solver. 

•  also  accurately  solves  time- dependent  diffusion  equations.  Note  that  the  treatment  of  diffusion  terms 
in  traditional  semi-discrete  DG  methods  are  non-trivial  because  of  the  so-called  “variational  crime”. 
In  DG-CVS,  the  inclusion  of  diffusion  terms  is  straightforward  and  simple  in  exactly  the  same  way 
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Figure  1:  Solution  elements  (SEs)  and  conservation  elements  (CEs)  in  the  x  —  t  domain.  Left:  solution 
elements  and  right:  conservation  elements. 


how  advection  terms  are  handled  by  simply  incorporating  the  diffusive  flux  into  the  space-time  flux 
[31].  No  extra  reconstruction  or  recovery  or  ad  hoc  penalty  and  coupling  terms  are  needed  to  ensure 
the  consistency  of  the  variational  form  for  diffusive  terms.  For  this  reason,  DG-CVS  is  conceptually 
simper  than  other  existing  DG  methods  for  diffusion  equations. 

The  DG-CVS  method  integrates  the  best  features  of  the  space-time  Conservation  Element /Solution  Element 
(CE/SE)  method  [34]  and  the  discontinuous  Galerkin  (DG)  method  [35].  The  core  idea  is  to  construct  a 
staggered  space-time  mesh  through  alternate  cell-centered  CEs  and  vertex-centered  CEs  (cf.  Fig.  1  (right)) 
within  each  time  step.  Inside  each  SE  (cf.  Fig.  1  (left)),  the  solution  is  approximated  using  high-order 
space-time  DG  basis  polynomials.  The  space-time  flux  conservation  is  enforced  inside  each  CE  using  the 
DG  discretization.  The  solution  is  updated  successively  at  the  cell  level  and  at  the  vertex  level  within  each 
physical  time  step.  For  this  reason  and  the  method’s  DG  ingredient,  the  method  was  named  as  the  space-time 
discontinuous  Galerkin  cell- vertex  scheme  (DG-CVS). 

DG-CVS  equally  works  on  higher  dimensions  on  arbitrary  grids.  Figure  2  shows  the  conservation  elements 
and  solution  elements  on  quadrilateral  meshes  and  triangular  meshes.  Obviously,  the  definitions  of  CEs  and 
SEs  on  higher  dimensions  are  analogous  to  that  for  1-D  meshes  (cf.  Fig.  1).  Figure  3  demonstrates  the 
resulting  dual  mesh  at  the  cell  level  and  the  vertex  level  for  both  rectangular  meshes  and  triangular  meshes, 
respectively. 

As  a  Riemann-solver  free  method,  DG-CVS  differs  from  other  Riemann- solver  free  central  schemes  [  6,  17] 
in  that  DG-CVS  employs  spacetime  discontinuous  Galerkin  polynomials  inside  the  solution  element  and 
therefore  is  arbitrarily  high  order  in  both  space  and  time.  Some  space-time  DG  methods  like  the  STE-DG 
(space-time  expansion  discontinuous  Galerkin)  scheme  [39]  also  approximates  solution  using  high-order  space- 
time  basis  polynomials  but  it  requires  some  numerical  flux  function  and  thus  is  subject  to  the  shortcomings 
associated  with  Riemann  solvers. 


2  Summary  of  Major  Accomplishments 

During  the  period  of  this  project,  we  have  accomplished  the  following: 

•  It  has  been  verified  that  the  DG-CVS  method  solves  3x3  MHD  model  equations  in  the  phase  space. 

•  The  method  has  been  successfully  extended  to  solve  the  ideal  ID  7-eq.  MHD  equations. 

•  The  method  has  been  successfully  extended  to  solve  the  ideal  2D  8-eq.  MHD  equations. 

•  It  has  been  verified  that  PowelTs  approach  by  adding  appropriate  source  terms  in  the  governing  equations 
is  effective  in  handling  the  V  •  B  =  0  condition  in  the  current  solution  framework. 
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cell  level  vertex  level 


Figure  2:  Conservation  elements  (CEs)  and  solution  elements  (SEs)  in  the  x  —  y  —  t  domain.  First  row:  CEs 
for  rectangular  meshes,  second  row:  CEs  for  triangular  meshes,  third  row:  SEs  for  rectangular  meshes,  and 
fourth  row:  SEs  for  triangular  meshes. 
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Figure  3:  Dual  meshes  for  the  solution  updating  at  the  cell  level  (in  red)  and  the  vertex  level  (in  black). 
Left:  rectangular  mesh  and  right:  triangular  mesh. 


3  Description  of  the  Spacetime  Discontinuous  Galerkin  Cell- Vertex 
Scheme 

To  illustrate  the  DG-CVS  formulation,  without  loss  of  generality,  we  consider  the  following  general  one¬ 
dimensional  time  dependent  partial  differential  equation  system 

+  dt=r  (Q\ 

dt  dx 

that  arises  from  some  physics.  Here,  u  is  the  solution  vector  containing  the  physical  quantities  to  be 
determined,  f  is  the  spatial  flux  vector.  Without  loss  of  generality,  a  source  term  r  is  also  included  in 
Eq.  (3).  Note  that  since  DG-CVS  solves  both  advection  and  diffusion  problems  in  an  identical  way,  f  may 
contain  both  advective  and  diffusive  fluxes.  Therefore,  f  can  be  functions  of  u  or  Vu  or  both.  Appropriate 
initial  and  boundary  conditions  must  also  be  provided  to  solve  Eq.  (3). 

3.1  Space-Time  Discontinuous  Galerkin  Formulation 

Following  the  idea  of  the  discontinuous  Galerkin  (DG)  method,  an  approximate  solution  uh  is  sought  within 
each  space-time  solution  element  (SE),  denoted  as  K.  When  restricted  to  the  SE,  uh  belongs  to  the  fi¬ 
nite  dimensional  space  U(K).  Assume  all  the  components  of  the  conservative  variables  inside  each  SE  are 
approximated  by  polynomials  of  the  same  degree,  i.e 

N 

(4) 

3  =  1 

where  {</q(x,  t)}^=1  are  some  type  of  space-time  polynomial  basis  functions  defined  within  the  solution 
element,  are  the  unknowns  to  be  determined  and  N  is  the  number  of  basis  functions  depending  on 

the  degree  of  the  polynomial  function.  Using  such  space-time  DG  solution  space,  the  solution  approximation 
can  be  arbitrarily  high  order  accurate  in  both  space  and  time. 

Following  the  Galerkin  orthogonality  principle,  multiply  (3)  with  each  of  the  basis  functions  fa  (i  = 

1,  2,  •  •  •  ,  N)  and  integrate  over  a  space-time  CE  to  obtain  the  weak  form 

r  /  r)uh  r)fh  \ 

where  Qk  is  the  space-time  conservation  element  (CE)  associated  with  the  solution  element  K.  Note  that 
the  conservation  element  is  identical  to  the  solution  element  except  for  the  volumeless  vertical  spike  in  the 
solution  element  as  shown  in  Fig.  ??.  The  space-time  flux  conservation  in  weak  form  as  in  (5)  is  for  each 
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individual  space-time  conservation  element.  Therefore,  the  current  method  can  be  considered  as  a  space- 
time  discontinuous  Galerkin  method.  As  can  be  seen,  the  space-time  flux  conservation  is  enforced  in  the 
variational  form  (5). 

Integrating  (5)  by  parts  results  in 


/ 

Jnj 


d(/)i 

dt 


nh 


dx 


d VL  =  j  4>iHh  ■  ndr 


(6) 


where 

=  ( ih ,  uh )  (7) 

is  the  space-time  flux  tensor  and  n  =  (nx,  nt)  is  the  outward  unit  normal  of  the  CE  boundary,  i.e.  T  = 
of  the  space-time  conservation  element  (CE)  under  consideration.  Note  that  the  partial  integration  is  also 
performed  on  the  time-dependent  term,  which  is  a  salient  difference  between  space-time  DG  methods  and 
semi-discrete  DG  methods.  As  can  be  seen,  the  formulation  in  (6)  contains  both  the  volume  integral  and 
the  surface  integral. 


3.2  Cell- Vertex  Solution  Updating  Strategy 

The  DG-CVS  inherits  the  core  idea  of  the  CE/SE  method  using  a  staggered  space-time  mesh  to  enforce 
the  space-time  flux  conservation.  However,  the  construction  of  the  staggered  space-time  slabs  in  DG-CVS 
deviates  from  the  CE/SE  method.  In  DG-CVS,  unknowns  are  stored  at  both  vertices  and  cell  centroids  of 
the  spatial  mesh,  and  the  solutions  at  vertices  and  cell  centroids  are  updated  at  different  time  levels  within 
each  time  step.  At  the  beginning  of  each  physical  time  step,  the  solution  is  assumed  known  at  the  vertices  of 
the  mesh,  either  given  as  the  initial  condition  or  obtained  from  the  previous  time  step.  Inside  each  new  time 
step,  the  solution  is  updated  in  two  successive  steps.  The  first  step  updates  the  solution  at  cell  centroids  at 
the  half-time  level  (£n+1/2)  based  on  the  known  vertex  solutions  at  the  previous  time  level  (£n).  The  second 
step  updates  the  solution  at  vertices  at  the  new  time  level  (£n+1)  based  on  the  known  cell  solutions  at  the 
previous  half-time  level  (£n+1/2).  The  same  process  is  repeated  for  new  time  steps. 

The  solution  updating  at  the  cell  level  or  the  vertex  level  is  based  on  the  key  equation  (6).  First  divide 
the  conservation  element  into  the  following  portions: 

•  the  interior  volume  Qk  where  the  solution  is  associated  with  the  space-time  node  at  the  new  time 
level. 

•  the  top  surface  T^0p  where  the  solution  is  also  associated  with  the  space-time  node  at  the  new  time 
level. 

•  the  side  surfaces  Tg^e  where  the  solution  is  associated  with  the  space-time  node  at  the  previous  time 
level. 

•  the  bottom  surfaces  T^0^  where  the  solution  is  also  associated  with  the  space-time  node  at  the  previous 
time  level. 

Correspondingly,  Eq.  (6)  can  be  rearranged  to  yield 


94>i 

dt 


H" • ndr ■ 


-L 

^rbott 


(f>iHh  •  ndr 


4>i&h  •  ndr 


**1,  2, 


(8) 


where  the  left  hand  side  contains  the  unknowns  and  the  right  hand  side  contains  the  known  information. 

Since  the  top  and  the  bottom  surface  of  the  CE  are  horizontal,  which  leads  to  •  n  =  uh  on  the  top 
face  and  •  n  =  -u^  on  the  bottom  face,  Eq.  (8)  can  be  simplified  to 


/ 

J  o  i 


d^uh  +  ^fh  +  4>irh 


<m- 


=  f  (j)iUh  ndr  -  [ 

"side  Jr 


/  0»u‘ 

Aop 


7;ii  dr  i,  2,  •  •  • ,  n 


iUhdT 


bott 


(9) 
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Figure  4:  Illustration  of  space-time  flux  conservation  on  a  1-D  cell-level  CE. 


To  further  illustrate  the  idea  of  enforcing  the  space-time  flux  conservation,  consider  the  conservation 
element  shown  in  Fig.  4.  Suppose  the  solution  at  the  spacetime  node  (m  +  \,n  +  |)  is  to  be  determined. 
Here,  m  and  n  represents  the  spatial  and  the  temporal  locations  of  the  space-time  node,  respectively.  The 
boundaries  of  the  CE  associated  with  the  spacetime  node  (m  +  \,n  +  is  divided  into  five  sections  T i, 
T 2 ,  r3,  r4  and  T5,  as  shown  in  Figure  4.  Among  these  sections,  T i  belongs  to  the  SE  associated  with  node 
(m  +  \,n  +  |)  whose  solutions  are  to  be  determined,  T2  and  T3  the  SE  associated  with  node  (m, n)  and 
r4  and  r5  the  SE  associated  with  node  (m  +  l,n).  The  interior  volume  of  the  conservation  element  is  also 
associated  with  node  (m  +  |,n  +  |).  Eq.  (9)  is  a  nonlinear  equation  system  for  each  space-time  node  which 
can  be  solved  by  the  Newton-Raphson  iterative  method. 

As  can  be  seen,  with  this  staggered  space-time  cell-vertex  solution  updating  strategy,  no  Riemann- solver- 
typed  flux  functions  are  needed  for  the  interface  flux.  There  is  no  “left  state”  and  “right  state”  when 
evaluating  inter-cell  fluxes  as  Riemann  solvers  do.  Our  Riemann- solver  free  DG-CVS  method  is  able  to 
capture  all  flow  features  (shock  waves,  contact  discontinuities,  etc.)  without  pathological  phenomena  such 
as  the  expansion  shock.  DG-CVS  equally  works  for  the  shallow  water  equations,  scalar  advect ion- diffusion 
equations  and  the  level  set  equation  in  addition  to  compressible  Euler  equations. 

3.3  Quadrature- Free  Implementation 

To  improve  efficiency,  the  integrals  on  the  left  hand  side  of  Eq.  (9)  are  evaluated  using  a  quadrature-free 
approach. 

To  allow  the  quadrature-free  implementation,  the  flux  and  source  vectors  must  be  expanded  in  terms  of 
the  basis  polynomials  similar  to  those  used  to  expand  the  conservative  variables,  namely, 

N  N  N 

fh  =  E  %h  =  E  sr fo' and  vh  =  E  s^j-  (10) 

j= 1  3  =  1  3  =  1 

where  N  is  the  number  of  basis  polynomials  of  one  degree  higher  than  those  to  expand  the  conservative 
variables  as  in  Eq.  (4).  The  requirement  of  using  basis  polynomials  of  degree  p  +  1  is  necessary  to  ensure 
the  accuracy  of  the  scheme. 

The  method  based  on  the  Cauchy- Kovalewski  (CK)  procedure [40,  41]  is  especially  suitable  for  our  pur¬ 
pose.  In  DG-CVS,  the  conservative  variables  are  expanded  using  the  Taylor  polynomials.  With  the  Taylor 
polynomials,  the  spatial  and  temporal  derivatives  of  the  solution  are  explicitly  available.  Using  the  CK 
procedure,  the  space-time  derivatives  of  the  flux  vectors  can  be  obtained  using  the  space-time  derivatives 
of  the  conservative  variables.  The  detailed  procedure  has  been  described  for  Euler  equation  in  our  earlier 
paper [33]  and  will  not  be  repeated  here. 


4  Major  Progress  #1:  Solving  the  MHD  Model  Equations 

The  present  Riemann- solver- free  method  has  been  first  applied  to  solve  the  following  3x3  MHD  model 
system  in  the  phase  space  which  exactly  preserves  the  MHD  singularities  [42]: 


ut  +  fx  =  duxx 


(ii) 
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pi  solution  of  Case  1  pi  solution  of  Case  1 


Figure  5:  Solutions  of  Case  1  by  solving  the  model  equations,  =  (— uo,  0,  0)  and  Ur  =  {uo,  2^0,  0)  with 

uo  =  1  at  t  =  0.0295085. 

with 

/  u  \  /  cu 2  +  u2  +  re2  \  /  /i  0  0  \ 

u  =  I  v  I  .  f  =  I  2  uv  I  ,  and  d  =  I  0  77  —  y  I 

\  w  )  \  2  uw  )  \  0  x  7]  ) 


where  the  longitudinal  viscosity  p  and  the  magnetic  resistivity  77  are  the  dissipative  coefficients,  and  the  Hall 
coefficient  y  is  a  dispersive  coefficient. 

The  above  quantities  are  related  to  the  MHD  system  via 


Bo, 


Bz 


a  =  I  —  —  1,  v  =  — -,  w  =  — ,  and  c  =  7  +  1 


H., 


where  a  =  \[ypjp  is  the  acoustic  wave  speed,  ca  =  ^p\Bx\  is  the  Alfven  wave  speed,  p  is  the  pressure,  p  is 
the  density  and  7  is  the  ratio  of  specific  heats. 

The  model  system  (11)  generates  three  wave  families  which  are  the  Alfven  wave,  the  slow  and  the  fast 
MHD  waves.  This  model  system  has  been  investigated  by  several  authors  [42,  43]  to  provide  insights  to 
improve  Godunov- typed  numerical  schemes  for  MHD  equations. 

The  two  test  cases  presented  in  this  section  are  taken  from  Ref.  [43]  where  conventional  and  modified 
Roe’s  scheme  were  used  to  solve  the  system.  In  all  the  cases,  c  =  3  and  the  computational  domain  [0, 1] 
is  evenly  divided  into  150  cells.  In  both  test  cases,  the  dissipative  term  on  the  right  hand  side  of  (11)  is 
assumed  to  be  zero. 


4.1  Case  1 

The  initial  discontinuity  is  a  fast  switch-off  expansive  shock  defined  as 


uL  =  {-Uo,  0,  0)  and  iir  =  {u0,  2 u0,  0) 

with  uo  =  1  in  this  case.  The  final  time  is  tm  0.0295085.  150  time  steps  are  run  to  reach  the  time  step. 

Figure  5  shows  u  vs.  x  and  v  vs.  x  solutions  from  the  pi  (second  order)  simulation.  As  can  be  seen, 
without  any  entropy  fix,  the  present  Riemann  solver  free  method  does  not  produce  expansion  shocks.  The 
comparison  between  the  analytical  solutions  and  the  present  solutions  demonstrate  the  accuracy  of  the 
present  method. 

4.2  Case  2 

The  initial  discontinuity  is  an  entropy- violating  shock  defined  as 
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pi  solution  of  Case  2 


pi  solution  of  Case  2 


x 


x 


Figure  6:  Solutions  of  Case  2  by  solving  the  model  equations.  ul  =  —  Ur  =  (—0.44,  2.4,  0)  at 
t  =  0.0376506024. 


U L  =  -UR  =  (-0.44,  2.4,  0) 

The  final  time  is  t  =  0.0376506024.  150  time  steps  have  been  run  to  reach  the  final  time. 

Figure  6  shows  u  vs.  x  and  v  vs.  x  solutions  from  the  second  order  simulation  along  with  the  analytical 
solutions.  Again,  the  agreement  is  very  well  and  the  overshoots/undershoots  around  the  discontinuity  region 
have  been  suppressed. 

5  Major  Progress  #2:  Solving  the  Ideal  ID  7-Eq  MHD  Equations 

Assuming  Bx  =  const.,  the  ideal  MHD  equations  reduce  to  a  7  x  7  system  where  the  state  vector  and  the 
flux  vector  are 


p 

pux 

pux 

pul+p+^Bf-Bl 

pUy 

puxuy  EXBy 

puz 

and  f  = 

puxuz  ExBz 

pE 

'LL X  (yPE  +  p  +  2  B  )  Bx  T  UyBy  T  uzBz) 

By 

LLx-^y  LLyBx 

_  B*  _ 

- 1 

H 

CQ 

K! 

53 

1 

CQ 

H 

53 

Since  Bx  is  constant,  the  divergence- free  condition  of  the  magnetic  field  is  automatically  satisfied. 

We  are  testing  the  test  cases  proposed  by  [5] .  The  problems  are  simulated  by  solving  the  7x7  1-D  Ideal 
MHD  Equations.  The  problem  to  be  solved  is  a  one-dimensional  Riemann  problem  on  the  computational 
domain  [—1, 1].  For  both  cases  below,  the  primitive  variables  are  defined  as  q  =  [p,u,v,w,p,  By,  BZ]T . 

5.1  Brio-Wu  Case  1: 

In  this  case,  the  initial  states  are  given  as 

qL  =  (1.000, 0, 0, 0, 1000, 1.0, 0.0)T  and  qR  =  (0.125, 0, 0, 0, 0.1,  -1.0, 0.0)T 

with  Bx  =  0  and  7  =  2. 

If  one  define 

p* =p+\  ibi2 

this  case  is  then  exactly  the  same  as  the  gas  dynamics  shock  tube  problem  except  the  equation  for  By. 
Figure  7  shows  the  numerical  solution  at  T  =  0.012  together  with  the  available  analytical  solution. 
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X 
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Figure  7:  Solutions  of  Brio-Wu  Case  1.  Bx  =  0.  T  =  0.012. 


It  can  be  seen  that  the  numerical  results  agrees  well  with  the  analytical  solution  except  for  the  slight 
overshoots/undershoots  around  discontinuities.  Note  that  no  any  type  of  limiting  procedure  is  employed  to 
obtain  the  current  solution  shown  in  Fig.  7. 

5.2  Brio-Wu  Case  2: 

In  the  second  case,  the  initial  states,  separated  at  x  =  0,  are  given  as 

qL  =  (1.000, 0, 0, 0, 1.0, 1.0, 0.0)T  and  qR  =  (0.125, 0, 0, 0, 0.1,  -1.0, 0.0)T 

with  Bx  =  0.75  and  7  =  2.  Figure  8  shows  the  unlimited  (left)  and  limited  (right)  solutions  at  T  =  0.2. 

As  can  be  seen,  the  complex  MHD  waves  including  the  compound  wave,  the  shock  wave,  contact  dis¬ 
continuity  and  rarefaction  waves  can  be  captured  by  the  present  Riemann-solver  free  method.  However,  the 
unlimited  solutions  exhibit  some  non-physical  overshoots  and  undershoots  near  discontinuities  (left  column 
in  Fig.  8).  These  overshoots/undershoots  can  be  suppressed  (right  column  in  Fig.  8)  using  a  minmod-typed 
limiting  procedure. 

6  Major  Progress  #3:  Solving  the  Ideal  2D  8-Eq  MHD  Equations 

The  8x8  ideal  MHD  equations  in  2-D  is  then  solved.  Here  the  state  vector  u  =  [p,  pux ,  puy ,  puz ,  pE ,  Bx ,  By ,  BZ]T 
and  the  flux  vectors  are 
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Figure  8:  Solutions  of  Brio-Wu  Case  2.  Bx  =  0.75.  T  =  0.2.  Left  column:  unlimited  solutions.  Right 
column:  limited  solutions. 
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The  2-D  case  solved  here  is  a  MHD  vortex  advection  problem.  The  solution  is  a  vortex  advected  with  the 
magnetic  field  diagonally  in  a  square  [—5,  5]  x  [—5,  5]  domain.  Periodic  boundary  conditions  are  assumed. 
The  computational  domain  is  discretized  by  a  32  x  32  rectangular  mesh.  The  solution  is  initialized  as  follows 

[19]: 


P 

u 

V 

w 


Bj; 

By 

Bz 


1.0 

l-0+^e°'5(1“ra)(— y) 

i.O  +  — e°'5(1-r2')x 
27 r 

0.0 

2n 

0.0 


where  r2  =  x2  +  y2. 

Figure  9  shows  the  magnetic  field  and  velocity  magnitude  field  at  T  =  0.0  and  T  =  10.0.  It  can  seen  that 
the  solution  fields  have  been  preserved  well  after  one  period’s  advection. 
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(a)  Bee-field.  Left  t  =  0.0,  right  t  =  10.0. 
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00 

(b)  B^-field.  Left  t  =  0.0,  right  t  =  10.0. 
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Ck 

(c)  Velocity  magnitude  field.  Left  t  =  0.0,  right  t  =  10.0. 

Figure  9:  MHD  vortex  advection  problem. 
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Randolph  Street  Suite  325,  Room  3112  Arlington,  VA  22203,  (703)  696-8429,  Fax  (703)  696-8450,  DSN 
426-8429,  fariba. f ahroo@af osr . af .mil. 

The  program  manager  has  been  changed  to  Dr.  Jean-Luc  Cambier  in  2015.  His  contact  email  is 

j  ean_luc  .  cambierQus  .  af  .  mil. 

11  Interactions/Transitions 

Conference  presentations  since  2013: 

1.  One  talk  at  the  SIAM  Conferences  on  Computational  Science  and  Engineering,  Salt  Lake  City,  Utah, 
March  2015.  (speaker:  S.  Tu). 

2.  One  talk  at  the  SIAM  Central  States  Conferences.  Rolla,  Missouri,  April  2015.  (speaker:  S.  Tu). 

3.  One  talk  at  the  AIAA  Orlando  SciTech  Conference  2015,  January  2015.  (speaker:  S.  Tu). 

4.  One  talk  at  the  AIAA  San  Diego  Summer  Conference  2013,  June  2013.  (speaker:  S.  Tu). 

5.  One  talk  at  the  SIAM  CSE  Conference,  Feb.,  2013,  Boston,  MA.  (speaker:  S.  Tu). 

12  Changes  in  Research  Objectives 


None. 


13  Change  in  AFOSR  Program  Manager 

The  program  manager  has  been  changed  to  Dr.  Jean-Luc  Cambier  in  2015. 

14  Extensions  Granted  or  Milestones  Slipped 

None. 

15  New  Discoveries 

None  patentable. 
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